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Q^' Abstract 

h>. ' In this work a simple method to enforce the positivity-preserving property for gen- 

o' 

Tij" ' eral high-order conservative schemes is proposed. The method keeps the original 

in 

scheme unchanged and detects critical numerical fluxes which may lead to neg- 

cn ■ 

^D ' ative density and pressure, and then imposes a cut-off flux limiter to satisfy a 

sufficient condition for preserving positivity. Though an extra time-step size condi- 
tion is required to maintain the formal order of accuracy, it is less restrictive than 



C^ , those in previous works. A number of numerical examples suggest that this method, 

when applied on an essentially non-oscillatory base scheme, can be used to prevent 
positivity failure when the flow involves vacuum or near vacuum and very strong 
discontinuities. 
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1 Introduction 

Compressible flow problems are usually solved by conservative schemes. High- 
order conservative schemes are suitable for simulating flows with both shock 
waves and rich flow features (acoustic waves, turbulence) since they are ca- 
pable of handling flow discontinuities and accurately resolve a broad range of 
length scales. One important issue of high-order conservative schemes is that 
non-physical negative density or pressure (failure of positivity) can lead to an 
ill-posed system, which may cause blow-ups of the numerical solution. While 
for some flrst-order schemes negative density or pressure can occur when a 
vacuum or near vacuum is reached, for higher-order conservative schemes pos- 
itivity failure can also occur due to interpolation errors at or near very strong 
discontinuities even though the flow physically is far away from vacuum. 



It is known that many flrst order Godunov-type schemes 2, ll6|, ISJ have the 

so called positivity-preserving property and can maintain positive density and 

pressure. It has been also proved that some second-order conservative schemes 
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5| are positivity-preserving with or without a more restrictive Courant- 
Friedrichs-Lewy (CFL) condition. For even higher-order conservative schemes, 
Perthame and Shu [llj] proved that, given a flrst-order positivity-preserving 
scheme, such as Godunov-type schemes, one can always build a higher-order 
positivity-preserving flnite volume scheme under the following constraints: (a) 
the cell-face values for the numerical flux calculation have positive density 
and pressure, (b) additional limits on the interpolation under a more restric- 
tive CFL-like condition. With a different interpretation of these constraints 
based on certain Gauss-Lobatto quadratures, positivity-preserving methods 
have been successfully developed for high-order discontinuous Galerkin (DG) 



methods 18| and weighted essentially non-oscillatory (WENO) finite volume 
and finite difference schemes [l9|, l20| . 



In this paper, we propose an alternative method to enforce the positivity- 
preserving property with a simple cut-off flux limiter. The flux limiter flrst 
detects critical numerical fluxes which may lead to negative density and pres- 
sure, then limits these fluxes to satisfy a sufficient condition for preserving pos- 



itivity. Unlike the approaches in 
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09 



20| , in which positivity-preserving and 



the maintenance of high order accuracy are considered simultaneously when 
designing the limiter, here we design the cut-off flux limiter to satisfy positiv- 
ity only, and then prove a posteriori the maintenance of high order accuracy 
under a time step restriction. It appears that, in our numerical experiments, 
a much less restrictive time-step size condition is sufficient for preserving pos- 
itivity without destroying overall accuracy. An advantage of the approach in 
this paper is that the cut-off limiter is directly applied to the numerical flux 
and it can be applied to arbitrary high-order conservative schemes. 



2 Method 

For presentation of the positivity-preserving flux limiters we assume that the 
fluid is inviscid and compressible, described by the one-dimensional Euler 
equations as 

where U = (p, m, E)'^, and F(U) = [m, pu^ + p,{E + p)u]'^. This set of equa- 
tions describes the conservation laws for mass density p, momentum density 



m = pu and total energy density E = pe + pM^/2, where e is the internal 
energy per unit mass. To close this set of equations, the ideal-gas equation 
of state p = (7 — l)pe with a constant 7 is used. Note that the density and 
pressure have the relations with the conservative variables as 

p(U)=p, p{V) = {^-l)(E-^-^y (2) 

It is easy to find that they are locally Lipschitz continuous, i.e. 

|p(U2)-p(Ui)|<L,||U2-Ui||, (3) 

|p(U2)-p(Ui)|<Lp||U2-Ui||, ifp(Ui)>0, p(U2)>0, (4) 

where Lp and Lp are Lipschitz constants. For 1 > ^ > 0, p(U) and p{U) have 
the properties 

p [(1 - ^)Ui + 0U2] = (1 - 0)p(Ui) + 0p(U2), (5) 

p [(1 - ^)Ui + 9V2] > (1 - %(Ui) + 9p{V2), if p(Ui), p(U2) > 0, (6) 

where Eq. i^ is straightforward and Eq. (jS]) is implied by the Jensen's in- 
equality since p(U) is a concave function. 

2.1 Finite-volume and finite- difference conservative schemes 

When Eq. ([T]) is discretized within the spatial domain such that Xj = zAx, i = 
0, ...,N, where Ax is the spatial step, a general explicit /cth-order conservative 
scheme with Euler-forward time integration can be written as 

U^i = Ur + A (f,_i/2 - F.+1/2) , (7) 

where the superscript n and n -\- 1 represent the old and new time steps, 
respectively, and A = At/Ax, where At is time-step size. Note that with the 



CFL condition 

CFL ■ Ax 

^t= , 

1 W -P Cjj^o;^ 



where c = J'yp/p is the sound speed and the CFL number < CFL < 1, one 
has the relation 



CFL 

^ = 7r-^ ^ — • 9 

[\U\ +C)max 



For a finite-volume scheme, U" and U"^^ are the cell averaged conservative 
variables on the cell i defined on the computational cell between {i — l/2)Aa; 
and {i + l/2)Ax, i.e. /j = [xi_i/2, Xj+1/2], Fi±i/2 = Fi±i/2 + 0(Ax''+^) are the 
numerical fluxes, which are based on the cell-face values Vi±i/2 reconstructed 
from the cell averages {Uj} and Fj±i/2 = F(Uj±i/2). 

For a finite-difference scheme, U" and U"^^ are the nodal values, and (Fj+1/2 — 
Fj_i/2)/Ax is a A;th order approximation to dF{\J)/dx at a; = Xi. Assume 
there exists a function H(a;) depending on Ax such that 



x+Ax72 

/ 

x-Ax/2 






then the same reconstruction procedure in a finite- volume scheme can be used 
to obtain the numerical fiuxes Fj±i/2 = iii±i/2 + 0{Ax''~^^) based on the 
cell- face values of H(x) reconstructed from its cell- average values F [\Jj] = 



/J-t/2' H(Oc?e/Ax. We refer to 



^Xj+Ax/2 

conservative finite difference schemes 



14j for the discussion of this formulation of 



2.2 Positivity preserving cut-off flux limiter 

The positivity-preserving property for the scheme Eq. ((Tj) refers to the property 
that the density and pressure are positive for U"^^ when U" has positive 
density and pressure. Since Eq. ([7]) can be rewritten as a convex combination 



Ur^ = ^ (ur + 2AF._V2) + ^ (Ur - 2AF.+V2 



2 

iu- + iu-, (11) 



a sufficient condition for preserving positivity is that Uj have positive density 
and pressure, i.e. g(Uf) > 0, where g represents p and p. Since the first-order 
Lax-Friedrichs fiux 



ir^F — 1 

' j+i/2 — 2 



F, + F,+, + {\u\+c)maA'U?-Ut 



+ly 



:i2) 



has the property g{\jf^'^) = g{\J'^ T'^^^ff^^^) > 0, under an additional CFL 
condition 

CFL < i (13) 



(see [18|), a straightforward way to ensure positivity is to hmit the magnitude 
of Fj_|_i/2 by utihzing the properties in Eqs. ([5]) and ([6]). The positive density 
is first enforced by: 

Cut-ofF flux limiter for positive density 



1. For all i: initialize 6'j^i/2 = 1) ^i+1/2 ~ ^- 



2. If p(U+) < e, , solve 9t,y, from (1 - 9t,y,)p{Vf^^) + 6ty2Pi^t) = ^p- 

3. If p(Ui;i) < tp, solve e-^^i^ from {l-ef^^i^)p{\J^^f)+ef^^i^p{\]^^-f) = e^. 

4. Set 9p^i+i/2 = min(6'+_^/2'^i+i/2)' ^i+1/2 = (l-6'p,i+i/2)Ff/i/2+^P,i+i/2Fi+i/2- 



Here, ep = min {10~^^, Pmj„}, where p^j„ is the minimum density in the initial 
condition, F*,,,^ ^ is the hmited flux, < di^i/2 ^ 1 are the hmiting factors 
corresponding to the two neighboring cells, which share the same flux Fj+i/2- 
After applying this flux limiter, Eq. flTTl) becomes 

= l^r + ^Ur . (14) 

Clearly, by Eq. (JS]), both U*'~ and U*' have positive density, so does U"^\ 
The positive pressure is further enforced by: 

Cut-off flux limiter for positive pressure 

1. For all i: initialize 6't';^i/2 ~ -'-' ^i+1/2 ~ -'-• 

2. If p(U:'+) < 6, , solve ^+ ,/2 from (1-0+ ,/,)p(Uf^'+) + e+ ,/2p(U*'+) = e^. 

3. Ifp(U*4ri) < ep, solve 0;;^/2 ^^^m {l-97^y^)p{Vf_^f) + eT-^^/^p{V*:^^) = ep. 

4. Set 9p,i+i/2 = min(6'^+ ^/2'^i+i/2)' ^^+1/2 = (l-^p,i+i/2)Ff/i/2+^P,i+i/2F*+i/2- 

Again, ep = min {10^^^,p|^j„}, where p^j„ is the minimum pressure in the 
initial condition, and F*^^ ,2 is the further limited flux. After applying this 
flux limiter, Eq. (TT4|) becomes 



ur^ = I (ur + 2AF*:i/2) + ^ (ur - 2Af*;i/ 

= l fur'' + ur'"") . (15) 



Clearly, by Eqs. (E]) and (jH]), both U**' and U**'+ have positive density and 
pressure, so does U"+^. Note that these limiters can be apphed at each sub- 



stage of a TVD Runge-Kutta 
Euler-forward time steps. 



13( 1 method, which is a convex combination of 



2.3 Consistency and accuracy 

Now we address two important issues for the cut-off flux limiter. First, tfie 
limited flux is a consistent flux since it is the convex combination of two 
consistent fluxes, i.e. the first-order Lax-Friedrichs flux U^ /g and the original 
high-order numerical flux F"^-,^^, which represents Fi+1/2 and F*_^^,2- Second, 
when the limiter is active, the difference between the original flux F°^^ /2 and 
the hmited flux F-![^ ,3 representing F*_^-^ ,2 and F*_|-,^ ,3 is 

l|F/^/2 - F°+i/2ll = (1 - ^s,i+i/2)||F°+i/2 - Fj_^i/2ll- (16) 

We only need to consider accuracy maintenance when 9gA+i/2 < 1, for other- 
wise the limiter does not take any effect. Without loss of generality we may 
assume 9g^i+i/2 = ^^i+i/2- ^^ ^^i^ situation we have g{\J°'^) < eg, in which 
U"'^ represents U^ and U*'^, and eg is negligibly small, and 

Since F°_^^ ,2 and F^ ,3 are both bounded in smooth regions, it is sufficient to 
show that the accuracy is not destroyed if the limiting factor satisfies 

l-^,,m/2 = 0(Ax'=+i), (18) 

a sufficient condition for which would be |5f(U°'^)| = 0{Ax''^^) and g{\Ji '"*") 
is bounded away from zero. 

Similar to Zhang and Shu 18], we assume the exact solution U(x) is smooth 
and g(Ui) > M, where Uj is either the cell-average (for the finite- volume 
scheme) or the nodal value (for the finite-different scheme) of the exact solution 



U(x) and M > is a constant. Since g{Vi) is obtained from a kth order 
approximation, one has g(Ui) > M — 0(Ax^^^) > AI/2 if Ax is sufficiently 
small, therefore 



5(U 



["•")= ^ 


/ 2\ -^ V 

{l-w)lJ, + w U,-— Ff/i/2 


/ 2A -^ ^ 


>(1 


7''m>o. 



(19) 



where 1 > u) > is a constant, under an extra CFL condition 



CFL < -. (20) 



Furthermore, one has 



ur=ur-2AF°+i/2 

= Uj ' + 2A (F^+1/2 - ^^+1/2 



Ur'+ + 2A F,^^,/2 - F.+1/2 + 0(Ax^+i), (21) 



where Fj+1/2 = Fj_|_i/2 for the finite-volume scheme, and Fj+1/2 = HJ4.1/2 for 
the finite-difference scheme. Let U| = Uj '^ + 2 A (F^ ^ — Fi+1/2) , and with 
Eqs. (|3]) and (jl]), one has 

l^(UD - ^(Ur)| < L,||Ur - Ulll = ©(Ax'^+i), (22) 

where Lg is the Lipschitz constant. Note that the first term of U| has positive 
density and pressure. For the second term, since the first-order Lax-Friedrichs 
fiux F^^/2 is a first order approximation to the exact fiux Fi^i/2, that is 
\\^i+i/2 - Fi+i/2|| = 0{Ax). With bounded ^(Uf^'+) from Eq. ([19]), one has 

Pm > ^^M - 0(Ax) > ^i^M > 



for sufficiently small Ax, according to Eq. ([3]), and furthermore p(U'l) > tp 
according to Eq. (j4]). Since g{\Jl) > eg and gCU"'^) < e^, i.e. p{\Jf) < tp while 
enforcing positive density and p(U*'^) < tp but p(U*''^) > tp while enforcing 
positive pressure, Eq. f l22|) leads to |f7(U"''^)| = 0{^x^^^). Hence, we have 
proved that the cut-off flux limiter preserves high-order accuracy. 

Note that, for given values of M and grid size, Eqs. f ITI?]) and f ET]) suggest 
that the errors introduced by the the cut-off flux limiter decrease with the 
time-step sizes. Also note that, the condition Eq. (120|) is less restrictive than 
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19l |20|, and is desirable for higher 



the time-step size conditions in Refs. 
computational efficiency. 

2.4 Assessment of accuracy 



As a simple way to test the accuracy of the present flux limiters, we consider 
the one- dimensional linear advection equation 

du du , , 

with initial condition u{x) > 0. Applying the cut-off flux limiter to preserve 
positivity results in the limiter (denoted as HAS) 



fUi/2 = ^(^r+1/2 - + <, = mill ; , 1 , 

M^in = min {< - 2AMr^i/2. <+i + 2AMr^i/2, IQ-^^} . (24) 

Here Mj+i /2 is the approximated upwind flux at the cell face i + 1/2. Note that 
only one of m" — 2\u^^-^^,2 or u^_^i + 2\u^^^,2 being negative will activate the 
limiter. The limiter of Zhang and Shu [l9|] (denoted as ZS) for Eq. (1231) can 



10 



be written as 



fUi/2 = ^K-+i/2 - <) +<, ^ = mill -^ , 1 , 






Umin = min { — ' /. ^^^^^, «+ 1/2' Wi+i/2' 10 ^H ' (25) 



where m^i/2 and ^i^i/2 are the high-order approximations of the cell-face 
values at m(xj_i/2) and 'u(xj+i/2) within cell i. For Eq. ([23]), one has m^^^ = 
Mj+i/2. Comparing Mmin in Eqs. fl24|) and fl25|) . it can be observed that the HAS 
limiter does not directly constrain the cell-face values to be non-negative. 

To further illustrate the accuracy of the HAS limiter and its relation to the ZS 
limiter, we compute the advection of a function u = 1 + 10~^ + cos(27rx) in do- 
main [0, 1] with a fifth-order conservative finite difference WENO-5 scheme 



t 



with third-order TVD Runge-Kutta time integration 13]. A periodic bound- 
ary condition is applied at x = and x = 1. The final time is t = 1, which 
corresponds to one period. This problem is computed on different grids with 
N = 50, 100, 200, 400 and 800 grid points. Figure [1^ shows the error dis- 
tributions for the results on 200 grid points. It can be observed that if the 
maximum admissible CFL number of 0.5 is used, the HAS limiter produces 
larger errors than the ZS limiter. However, the HAS limiter is already as accu- 
rate as the ZS limiter when a smaller CFL number of 1/12, which corresponds 
to the maximum admissible value for the latter, is used. If the time-step size 
is decreased further, errors produced by the ZS limiter do not change con- 
siderably, whereas the errors produced by the HAS limiter decrease further. 
This behavior is also shown in Fig. [Jo for the evolution of the L^q error with 
decreasing grid size. Here, the time-step size At = 0.5Ax^/^ is used to keep the 
spatial errors dominant. Note that Fig. [T)d clearly shows that the theoretical 

11 



order of accuracy is achieved. 

2.5 Extension to multiple dimensions 

To present the extension of the positivity-preserving flux hmiters to multiple 
dimensions we consider the two-dimensional Euler equation 

m ^ gF(U) ^ dG{V) ^ ^ 
dt dx dy 

where U = (p, pw, pv, E)'^, F(U) = [pu, pu^ + p, puv, {E + p)u]'^ and G(U) = 
[pv, pvu, pv'^+p, {E+p)v]'^. Compared to the one- dimensional equation Eq. ([1]), 
the momentum density is pv = {pu, pv), where u and v are velocities in the x 
and y directions, respectively, and the total energy density is -E = pe + pv^/2. 
As an extension of Eq. ( fTTj) . the conservative scheme for Eq. (!26|) can be 
rewritten as a convex combination 



+ ^ (Ui:,. + 2A,G,,,_i/2) + f (Ui:,- - 2A,G,,+i/2) , (27) 



where A^ = At/Axax and Ay = At/ Ay ay, ax + ay = 1, with a^ > and 
ay > being partitions of the contribution in the x and y directions. A simple 
way to obtain this partition is to set ax = ay = 1/2 as in Zhang and Shu 



m 



20|. Another straightforward way to determine ax and ay = 1 — ax is 

Tx (|M|+C)max (|^|+c)max /^t 

O^x = ■ , Tx = -r , Ty = . {2t 

Tx + Ty Ax "^ Ay 



Note that, since the time-step size for integrating Eq. fl27p is given by 

At = ^. (29) 
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one has the relation 

CFL CFL 

(|M|+C)max (|f|+c)max 

which gives an extended form from Eq. (Q. Also note that, since the compo- 
nents in Eq. ( 127|) and Eq. ( ITT]) have the same form, it is straightforward to 
implement the positivity-preserving flux limiters in a dimension-by- dimension 
fashion. 



3 Test cases 

In the following, we illustrate that a number of typical numerical test cases, 
where the original high-order conservative schemes fail, can be simulated by 
using the proposed positivity-preserving flux limiters. For the first type of 
cases involving vacuum or near vacuurn, the flux limiters are combined with 



the finite difference WENO-5 scheme 



7|], which is a shock-capturing scheme 



with fifth-order accuracy for smooth solutions. For the second type of cases 
involving very strong discontinuities, the flux limiters are combined with the 
WEN0-CU6-M1 scheme ^, which can be used for implicit large eddy simula- 
tion (LES) of turbulent flow and has sixth-order accuracy for smooth solutions. 
For both variants of the WENO schemes the Roe approximation is used for the 
characteristic decomposition at the cell faces, the Lax-Friedrichs formulation is 
used for the numerical fluxes, and the third-order TVD Runge-Kutta scheme 
is used for time integration [l3|. If not mentioned otherwise, the computations 
are carried out with a CFL number of 0.5. 
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3.1 One- dimensional problems involving vacuum or near vacuum 

Here we show that the proposed method passes two one-dimensional test prob- 
lems involving vacuum or near vacuum: the double rarefaction problem [5|, 

^1, 



where a vacuum occurs, and the planar Sedov blast-wave problem 

where a point-blast wave propagates. For the first problem, the initial condi 

tion is 

_ ;i,-2,0.1) if <a; <0.5 



1,2,0.1) ifl>x>0.5 



-) "1 



Ax = 2.5 X 10 ^ and the final time is t = 0.1. For the second problem, the 
initial condition is 



_ (1,0,4 X 10-13) if <x <2-0.5Ax, 2 + 0.5Ax < x < 4 
(p>M,p) \ 

(1, 0, 2.56 X 10^) if 2 - 0.5Ax < x < 2 + 0.5Ax 



Ax = 5 X 10 3 and the final time is t = 10 3. 

Figure |5]gives the computed pressure, density and velocity distributions, which 
show good agreement with the exact solutions. Although a vacuum occurs in 
the solution of the double rarefaction problem, the results still exhibit accu- 
rate density and pressure profiles in the rarefaction-wave regions. As a vacuum 
occurs, the solution at the center of the domain strictly speaking has no phys- 



ical meaning. Note that compared to Zhang and Shu 20| (see their Fig. 5.1 
(right)) for the planar Sedov blast- wave problem a slightly sharper blast wave 
is obtained in the present results. This may be due to the fact that Zhang 



and Shu 



20| have modified the original Lax-Friedrichs flux to use a single 



maximum signal speed other than the respective maximum eigenvalues. 
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3.2 Two-dimensional problems involving vacuum or near vacuum 

We consider two two-dimensional problems involving vacuum or near vacuum. 
The first problem is the two-dimensional Sedov problem which has been stud- 

nn 

ied in Zhang and Shu [18|, 120| . 



The computation is performed on the domain 
[0,0] X [1.1, 1.1], where a high pressure region occupies the computation cell 
at the lower-left corner. The initial condition is given by 

1,0,0,4x10-^3) ifx>Ax, y>Ay 



{p,u,v,p) 



:i'0'0'^xlO^) else 



where Ax = Ay = 1.1/160. The final time is t = 1.0 x 10"^. A reflective 
boundary condition is applied at the lower and left boundaries, and an outflow 
condition is applied at the right and upper boundaries. Figure |3] gives the 
computed density profiles. One can observe that these results are in very good 
agreement with the exact solution. 

The second problem is the Mach-2000 jet problem, which has been computed 



in Zhang and Shu 18 



19 



201 ]. The computation is performed on the domain 
[0,1] X [0,0.25], Initially, the entire domain is filled with ambient gas with 
{p,u,v,p) = (0.5,0,0,0.4127). A reflective condition is applied at the lower 
boundary, an outflow condition is applied at the right and upper boundaries, 
and an inflow condition is applied at the left boundary with states (p, u, v, p) = 
(5,800,0,0.4127) iiy < 0.05 and {p,u,v,p) = (0.5,0,0,0.4127) otherwise. A 
CFL number of 0.25 is used and the flnal time is 0.001. Since 7 = 5/3 is 
used, the speed of the jet 800 gives about Mach 2100 with respect to the 
sound speed in the jet gas. Figure H] gives the computed density and pressure 



profiles in logarithmic scale. One can observe 



;hat these results are in very 



good agreement with those in Zhang and Shu 18| (their Fig. 4.6) computed 
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with the same resolution. 

3.3 One- dimensional problems involving very strong discontinuities 

We show that, combined with the proposed flux limiters, the WEN0-CU6-M1 
scheme passes two one-dimensional test problems, which cannot be computed 
with the original scheme without limiting, involving very strong discontinu- 
ities: the two blast-wave interaction problem 



17j , and the Le Blanc problem 



{P,u,p) 



10l |. The latter is an extreme shock-tube problem. For the first problem, the 
initial condition is 

f (1,0, 1000) ifO<x<0.1 

1,0,0.01) if0.1<x<0.9, 
[(1,0,100) ifl>x>0.9 
Ax = 2.5x 10"'^, and the final time is t = 0.038. Reflective boundary conditions 
are applied at both x = and x = 1. The reference "exact" solution is a high- 
resolution numerical solution on 3200 grid points calculated by the WENO- 



CU6 scheme 



6|]. For the second problem, the initial condition is 
1,0, 1 X 10-1) if0<x<3 



{p,u,p) = 
7 = 5/3, Ax = 9/800 and the final time is t = 6. 



10-3, 0, f X 10-1°) if3<x<9 



Figure [5] gives the computed pressure, density and velocity distributions, al- 
though at relatively low resolution, which show a good agreement with the 
exact or reference solutions. The magnitudes of the small over-shoots (see Fig. 
|5]^left)) and the small errors at the shock position (see Fig. [5]^ right)) decrease 
when the resolution is increased (not shown here). For the two blast-wave 
interaction problem the present results are comparable to those obtained by 
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the WENO-CU6-M2 scheme ^ at the same resolution (see their Fig. 3). Note 
that the WENO-CU6-M2 scheme stabihzes for very strong discontinuities in 
a different way, but still cannot compute the Le Blanc problem. 

3.4 Two-dimensional problems involving very strong discontinuities 



We first consider the problem from Woodward and Colella 17| on the double 

Mach reflection of a strong shock. A Mach 10 shock in air is reflected from 

the wall with incidence angle of 60°. The initial condition is 

r (1.4, 0,0,1) if ?/< 1.732(a;- 0.1667) 

{p,u,v,p) = \ 

[ (8, 7.145, -4.125, 116.8333) else 

and the final time is t = 0.2. The computational domain for this problem is 

[0,0] X [4,1]. Initially, the shock extends from the point x = 0.1667 at the 

bottom to the top of the computational domain. Along the bottom boundary, 

at y = 0, from x = to x = 0.1667 the post-shock conditions are imposed, 

whereas a reflective condition is set from x = 0.1667 to x = 4. Inflow and 

outflow conditions are applied at the left and right boundaries, respectively. 

The states at the top boundary are set to describe the exact motion of a Mach 

10 shock. Figure [6] shows the pressure and density contours of the solution on 

a 240 X 60 grid. Note that compared to the results obtained by WEN0-CU6- 



M2 



J] (their Fig. 4) a good agreement is observed. Especially, both predict 



a strong near-wall jet, which is usually smeared in the previous computations 



with the same resolution 
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l8|,l6|. 



We then consider a shock-bubble interaction problem, when a Mach 6 shock 
wave in air impacts on a cylindrical helium bubble. Air and helium are treated 
as the same ideal gas fluid for simplicity. Numerical computations for this 
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problem can be found in Bagabir and Drikakis l|. The initial conditions are 

' {p = l,u = —3, V = 0,p = 1) pre-shocked air 

{p = 5.268, u = 2.752, v = 0,p = 41.83) post-shocked air , (31) 

^ (p = 0.138,M = — 3,f = 0,p = 1) helium bubble 

and the final time is t = 0.15. The computational domain for this problem is 
[0,0] X [1,0.5]. Initially, the shock wave is at x = 0.05, and the half helium 
bubble of radius 0.15 is at (0,0.25). Note that a frame velocity u = —3 is 
applied to keep the bubble approximately in the center of the computational 
domain. Refiective conditions are applied at the lower and upper boundaries, 
an outfiow condition is applied at the right boundary, and an infiow condition 
is applied to the left boundary with the post-shocked state. Figure [7] shows 
the pressure and density contours of the solution on a 200 x 100 grid. These 
results show a fairly good agreement with those in Bagabir and Drikakis 1| 
(their Fig. 6) at the same resolution. The secondary refiected shock wave 
and triple-wave configurations are calculated with good resolution. Note that 
since the WEN0-CU6-M1 scheme has smaller numerical dissipation than the 
MUSCL scheme used in Bagabir and Drikakis [1], the present results show a 
less smeared bubble interface and more detailed structures near the triple- wave 
region. 

4 Concluding remarks 

In this paper we have proposed a very simple method to enforce the positivity- 
preserving property for general high-order conservative schemes. The method 
first detects critical numerical fiuxes which may lead to negative density and 
pressure, then limits the fiuxes to satisfy a sufficient condition for preserving 



positivity. Though an extra time-step size condition is required to maintain 
the formal order of accuracy, it is less restrictive than those in previous works. 



In addition, since the method uses the general 



similarly as the approaches of Zhang and Shu 20| , it can be applied to flows 



brm of a conservative scheme. 



with a general equation of state and source terms in a straightforward way. 
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Figure 1. Linear advection problem at t = 1: (a) Error distribution vs. time-step 
sizes on 200 grid points; (b) Evolution of the L^o error with decreasing grid size.. 
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Figure 2. One-dimensional problems involving vacuum or near vacuum: (left) double 
rarefaction problem; (right) planar Sedov blast-wave problem. 
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Figure 3. Two-dimensional Sedov problem: (left) 10 density contours from to 6; 
(right) density profile along y = 0. 




Figure 4. Mach-2000 jet problem: (upper) 30 density contours of logarithmic scale 
from -4 to 4; (lower) 30 pressure contours of logarithmic scale from -1 to 13. 
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Figure 5. One-dimensional problems involving very strong discontinuities: (left) two 
blast wave problem; (right) Le Blanc shock-tube problem. 
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Figure 6. Double-Mach reflection of a Macli 10 sliock wave at t = 0.2: (upper) 30 
pressure contours from 0.92 to 520; (lower) 30 density contours from 1.73 to 21. 
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Figure 7. Shock-bubble interaction problem at t = 0.15: (left) 30 pressure from 0.9 
to 62; (right) 30 density contours from 0.5 to 8. 
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